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ASSESSMENT OF STATIC DELAMINATION PROPAGATION CAPABILITIES IN 
COMMERCIAL FINITE ELEMENT CODES USING BENCHMARK ANALYSES 


1 2 

Adrian C. Orifici , Ronald Krueger 

ABSTRACT 

With the increasing implementation into commercial finite element (FE) codes of 
capabilities for simulating delamination growth in composite materials, the need for 
benchmarking and assessing these capabilities is critical. In this study, benchmark analyses 
were performed to assess the delamination propagation simulation capabilities of the VCCT 
implementations in Marc 2008rl and MD Nastran R3 v2008.0 in solution sequences 
SOL 400 and SOL 600. Benchmark delamination growth results for Double Cantilever 
Beam (DCB), Single Leg Bending (SLB) and End Notched Flexure (ENF) specimens were 
generated using a numerical approach. This numerical approach was developed previously, 
and involves comparing results from a series of analyses at different delamination lengths to 
a single analysis with automatic crack propagation. Experimental and analytical benchmark 
delamination growth results were also taken from the literature for a second DCB specimen. 
Specimens were analyzed with three-dimensional and two-dimensional models, and 
compared with previous analyses using Abaqus® with the VCCT implemented. The results 
demonstrated that the VCCT implementation in Marc and MD Nastran was capable of 
accurately replicating the benchmark delamination growth results. The analyses in Marc 
and MD Nastran™ were significantly more computationally efficient than previous analyses 
conducted in Abaqus®. This was due to a lack of convergence issues, and a solution process 
that maintained the use of large time increments. The Marc™ and MD Nastran™ solvers 
encountered problems for some models with determining the appropriate crack growth 
direction, which required overriding the default automatic procedure. The results also 
demonstrated that the use of the numerical benchmarks was efficient and offers advantages 
over benchmarking using experimental and analytical results. 


1. INTRODUCTION 

One of the most common failure modes for composite structures is delamination [1-4]. To 
characterize the onset and propagation of delamination, the use of fracture mechanics has become 
common practice over the past two decades [1, 5, 6]. The strain energy release rate, Gr, is typically 
used as a measure of the driving force for delamination growth in composite laminates. Depending 
upon external loading conditions, Gt can be any combination of the three strain energy release rate 
components, G/, G//, and G///, illustrated in Figure 1. The mode I component, G/, arises from loading 
acting normal to the delamination plane, Gn arises from in-plane shear stresses acting perpendicular 
to the delamination front, and Gm arises from in-plane shear stresses acting parallel to the 
delamination front. To predict delamination onset or propagation, the total strain energy release rate 
is compared to the interlaminar fracture toughness, G c , which is dependent on the relative 
proportions of the mode components. Due to the availability of test methods for characterizing 
mode I, II and mixed mode I/II delamination, efforts [7-9] have focused on evaluating the 
dependence of G c on this range of mode mix. Such a quasi-static mixed-mode I/II fracture criterion 


1 RMIT University, GPO Box 2476, Melbourne, Victoria 3001, Australia. 

^ National Institute of Aerospace (NIA), 100 Exploration Way, Hampton, VA 23666, resident at Durability, Damage 
Tolerance and Reliability Branch, NASA Langley Research Center, MS 188E, Hampton, VA, 23681. 
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is determined by plotting G c versus the mixed-mode ratio, Gn/Gr. The mixed-mode ratio is 
determined from data generated using pure mode I Double Cantilever Beam (DCB) (Gn/Gi=0), 
pure mode II End-Notched Flexure (ENF) (G//G 7 =l), and Mixed-Mode Bending (MMB) tests of 
varying ratios. Examples are shown in Figure 2a for T300/914C and Figure 2b for C12K/R6376 
[10, 1 1]. A curve fit of the data is performed to determine a mathematical relationship between G c 
and Gn/Gr. [12, 13]. An interaction criterion incorporating mode III was recently proposed by 
Reeder [14], Application of Reeder’s 3D criterion is restricted due to the absence of a reliable 
means for measuring Gm c . Although, a standardized test for measuring this property is currently 
being investigated [15, 16], and several candidate specimens have been proposed in the literature 
[17,18]. 

The virtual crack closure technique (VCCT) is widely used for computing energy release 
rates, based on results from continuum (2D) and solid (3D) finite element (FE) analyses, and to 
supply the mode separation required when using mixed-mode fracture criteria [19, 20]. The VCCT 
has been used mainly by scientists in universities, research institutions and government laboratories, 
and is usually implemented in their own specialized codes or used in post-processing routines in 
conjunction with general purpose FE codes. An increased interest in fracture mechanics to assess 
the damage tolerance of composite structures in design and certification has also renewed interest in 
the VCCT. As such, the VCCT has been implemented into the commercial FE codes Abaqus® 1 [21] 
as well as MD Nastran™ [22] and Marc™ 2 [23], among others. The implementation into the 
commercial FE code SAMCEF™ 3 [24] is a mix of the VCCT and the Virtual Crack Extension 
Method suggested by Parks [25]. 

As new approaches for analyzing composite delamination are incorporated into FE codes, 
the need for comparison and benchmarking becomes important. The implementation of 
technologies into FE codes involves numerical parameters, which can be unique to each code, and 
need to be understood and calibrated for any analysis. Experimental results used as benchmarks are 
valuable, but are complicated by aspects such as fiber bridging, crack branching, and experimental 
variance. Analytical results are also useful as benchmarks, but are not available for all specimen 
types and configurations, and can become complicated when dealing with aspects such as 
nonlinearity and three-dimensional (3D) effects. In response, a numerical benchmarking approach 
was developed in which results for delamination growth are first generated from a series of static 
analyses with different delamination lengths [26]. These benchmark results are then compared to 
simulation of delamination propagation using a single analysis. In previous work, this numerical 
benchmarking approach was applied to DCB and Single Leg Bending (SLB) specimens to assess 
the implementation of VCCT in Abaqus®/Standard [26]. 

In this study, the delamination propagation simulation capabilities of the commercial FE 
codes Marc™ 2008rl and MD Nastran™ v2008.0 with the VCCT were assessed. Benchmark 
delamination growth results for DCB, SLB and ENF specimens were generated using a previously 
developed numerical approach [26]. Experimental and analytical benchmark delamination growth 
results were taken from the literature for a second DCB specimen [27]. Specimens were analyzed 
using two-dimensional (2D) plane strain and full 3D models and results were compared to previous 
analyses using Abaqus®/Standard with the VCCT implemented. 


1 Abaqus® is manufactured by Dassault Systemes Simulia Corp. (DSS), Providence, RI, USA 

2 MD Nastran™, MSC Nastran™, MSC Patran™, Marc™ and Mentat™ are manufactured by MSC.Software Corp., 
Santa Ana, CA, USA. NASTRAN® is a registered trademark of NASA. 

3 SAMCEF™ is manufactured by Samtech, Liege, Belgium. 
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2. SPECIMEN DESCRIPTION 

For the current investigation, specimens were selected to investigate single mode (I and II) 
and mixed-mode I/II delamination growth with unidirectional and multi-directional laminates, using 
benchmark results from numerical and experimental studies. Two DCB specimens (labeled “DCB- 
1” and “DCB-2”), one SLB specimen and one ENF specimen were chosen, as shown in Figure 3, 
Figure 4 and Figure 5, respectively. The specifications for the DCB-1 and SLB specimens 
originated in previous work, in which experimental and numerical studies were performed and the 
critical strain energy release rates evaluated [28-31]. In this work, the DCB-1 and SLB specimens 
were selected as they were used to generate numerical benchmark delamination growth results 
previously [26]. The ENF specimen is the three-point bending 3ENF variant, and the specifications 
were set based on commonality with the DCB-1 specimen. The DCB-2 specimen (shown in 
Figure 3) was selected as benchmark experimental and analytical results were available in the 
literature [27]. The DCB and ENF specimens are used to determine the mode I and II interlaminar 
fracture toughness, G/c [7] and Gnc [9], respectively. The SLB specimen is used to evaluate the 
mixed-mode I/II interlaminar fracture toughness [30, 32], 

The DCB-1 and ENF specimens consisted of T3 00/ 1076 graphite/epoxy with a 
unidirectional layup and the DCB-2 specimen consisted of T800/924 graphite/epoxy with a 
unidirectional layup. The SLB specimen used C12K/R6376 graphite/epoxy with a multi-directional 
layup of stacking sequence [±3 0/0/- 3 0/0/3 O/O4/3 0/0/-3 0/0/-3 0/3 0//- 3 0/3 0/0/3 0/0/ -3 O/O4/- 

30/0/30/0/±30], where the double slash (//) denotes the location of the delamination. The material 
properties are given in Table 1. In general, mode I, mode II, and mixed-mode tests are performed on 
unidirectionally reinforced laminates where all fibers are aligned with the 0° direction. This means 
that delamination propagation occurs at a [0/0] interface and delamination propagation is parallel to 
the fibers. Although this unidirectional layup is desired for standard test methods to generate 
fracture toughness data, delamination propagation between layers of the same orientation will rarely 
occur in real structures. Therefore, effects such as fiber bridging that are observed in DCB tests 
may not be be present during mode I-dominated delamination growth between plies of dissimilar 
orientation. 


3. METHODOLOGY 
3.1 Fracture Criteria 

Linear elastic fracture mechanics analysis of delamination in composite laminates involves 
determining the total strain energy release rate, Gt, and the individual components G/, Gn and Gm, 
shown in Figure 1 . The onset of delamination growth is predicted using the failure index: 


where Gt is the sum of all mode components and G c is the interlaminar fracture toughness. The 
fracture toughness is dependent on the relative proportions of the mode components, or the mode 
mix. A 2D relationship between G c and modes I and II was suggested by Benzeggah and Kenane 
[13] and is given by: 


G c = G Ic + { G llc 



( 2 ) 
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where G/ c and Gn c are determined experimentally from DCB tests [7] and ENF tests [9], 
respectively, and the 2D mixed-mode exponent r] is determined from MMB tests of varying ratios 
of Gi and Gu [8]. Figure 2a and 2b show the determination of these parameters from experimental 
data for the materials applied in this work, where r] was determined by a curve fit using the 
Levenberg-Marquardt algorithm in KaleidaGraph™ 3 graphing and data analysis software [33], 

Recently, Reeder [14] suggested a modification to Equation 2 that incorporates the mode III 
component, which is given by: 


G c = G Ic + { G IIc 




( 3 ) 


As no standards currently exist for determining Gm c , or a mixed-mode exponent involving mode III, 
in this work Gm c was taken as G// c and rj was taken as the value characterized from 2D mode I/II 
data as shown in Figure 2. 

3.2 Virtual Crack Closure Technique 
3.2.1 Theory 

The VCCT [19, 20] is based on the assumption that the energy released in extending a crack 
by a small amount, A a, is equivalent to the work necessary to close the crack to its original length. 
In an FE analysis using the VCCT, the three strain energy release rate components are calculated at 
a crack front node by 


G = 


F-6 
2A ’ 


( 4 ) 


where G is the column vector of strain energy release rate components, F consists of the forces at 
nodes along the delamination front, 6 is a row vector that consists of the relative displacements of 
the node pairs behind each corresponding crack front node, and A is the surface area created by 
crack growth. Equation 4 requires the calculation of a local crack front coordinate system and 
modification to account for arbitrary element sizes. The VCCT is applicable for 2D or 3D analysis 
with linear and quadratic elements [20]. 

Previous investigations have shown that the VCCT results can be mesh-dependent for 
delaminations between two dissimilar solids, such as between plies of different orientation. This is 
due to the oscillatory nature of stresses in the vicinity of a crack front at a bi-material interface. To 
avoid mesh dependency issues, it is recommended to use consistent crack tip element lengths 
between analyses that are being compared. Further detail on this issue is provided in Reference [20]. 
In the current work, the finite element meshes were all based on those used in previous analyses 
[26], where mesh convergence studies had been conducted. 


3 KaleidaGraph™ is manufactured by Synergy Software, Reading, PA, USA. 
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3.2.2 VCCT in Marc™ 

The VCCT is implemented into Marc as a procedure to determine the strain energy release 
rate distribution at a crack front [23]. This is done at a single node in a 2D model, or along a series 
of nodes in a 3D model. Multiple cracks can also be modeled. The user only needs to define the 
crack front nodes, and the solver determines the appropriate nodes, forces and areas to use for the 
crack growth calculation. This implementation follows the description given by Krueger [20], which 
accounts for a crack front of arbitrary shape. 

For propagating a crack front, there are several options: growth by remeshing, where the 
local mesh around the crack tip is moved and the surrounding mesh is adapted; growth along 
element edges, where the elements at the crack front are disconnected by creating duplicate nodes 
and modifying the element connectivity; growth by releasing constraint, where a fixed contact or tie 
restraint is released. In the current work, only the latter approach was applied, as this was the only 
approach that could be used with 3D solid elements. 

For crack growth by releasing constraints, the user defines two surfaces that are initially 
bonded with a “glued” [23] constraint, or a set of multi-point constraints (MPCs) between nodes 
across an interface. In the current work, the contact-based approach was used, which is similar to 
the approach implemented within Abaqus® [21]. It is possible to define a completely bonded 
interface, and then to select certain nodes that are unbonded, or have the bond function 
“deactivated” [23], so that a pre-crack can be modeled. An example of this is shown in Figure 6, 
where a specimen with a pre-crack (shown in Figure 6a) is modeled using only a bonded interface 
(shown in Figure 6b), or with a combination of a bonded and unbonded contact interface (shown in 
Figure 6c). 

At the end of every nonlinear analysis increment, the strain energy release rates are 
calculated using the VCCT. Crack growth onset is detected using Equation 1, with either single- 
mode or mixed-mode criteria for G c . For single-mode criteria, three equations are used to consider 
each mode independently, and for mixed-mode, the Reeder (Equation 3) or Power law criteria are 
available [14]. The crack is “grown” by releasing the constraint at the crack front node, which 
extends the crack to the next node in the bonded contact region. If crack growth is detected at the 
end of the increment, the appropriate nodes are released, and the increment is restarted. Restarting 
the increment is repeated until Equation 1 is no longer satisfied at any crack front node, and is a 
critical step in ensuring that enough crack growth occurs within an increment This aspect is detailed 
further in the Discussion section. 

Another aspect that can become important is the way in which the VCCT algorithm locates 
the most appropriate node in the intact region to “grow” the crack, once propagation has been 
detected at a crack front node. In Marc , this is based on the definition of the crack growth 
direction, where once this direction is determined, the solver will select the closest appropriate node 
from those available in the intact region. The default option in Marc is an automatic determination 
of the crack growth direction based on a maximum principal stress criterion, which defines the 
crack growth as being normal to the direction of greatest tension [23]. There are other options for 
specifying the crack growth direction that the user can select instead of the default, where the 
growth direction is instead aligned with either the most critical mode, with the mode I component, 
or with a fixed user-defined vector. The effect of this option is discussed further in the results for 
each specimen, and further detailed in the Discussion section. 
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3.2.3 VCCT in MD Nastran™ SOL 400 

The VCCT is implemented into MD Nastran™ solution sequence SOL 400, and follows the 
same general approach as outlined above for Marc . The SOL 400 solution sequence is an 
integration of the traditional nonlinear solver (SOL 106) with advanced nonlinear technologies [23]. 
This includes functionality developed originally in Marc , such as contact and the VCCT 
technology for crack growth modeling. Critically though, not all aspects of the Marc 
implementation within a given software version are available in the MD Nastran version of the 
same release date [23, 22, 34], As a key example of this limitation, in the SOL 400 version applied 
in this work (version 2008.0), the VCCT technology did not incorporate mixed-mode crack growth. 
Another important example is that only the default option for determining the crack growth 
direction could be specified as input. 

3.2.4 VCCT in MD Nastran™ SOL 600 

The VCCT is available in MD Nastran™ solution sequence SOL 600. The SOL 600 solution 
sequence is a “wrapper” application for Marc , in which the user generates an MD Nastran input 
file, and the solver converts this into a Marc input file, runs the analysis natively in Marc , and 
translates the Marc output into MD Nastran format [22], However, despite running Marc 
natively, this solver can be limited by the translation process, which includes options available in 
Marc not being easily available within the MD Nastran input format, and some difficulties in 
“mapping” functionality between the two codes. As with the SOL 400 solver, a key example of the 
limitations of the SOL 600 solver in the version applied in this work was that the VCCT technology 
did not incorporate mixed-mode crack growth, and could only be used with the default option for 
determining the crack growth direction. While the SOL 400 solution sequence is exclusive to the 
MD Nastran analysis package, the SOL 600 solution sequence is also available within MSC 
Nastran . As such, references made in this work to MD Nastran SOL 600 are identically 
applicable to MSC Nastran™ SOL 600. 

4. FINITE ELEMENT MODELING 

Typical 2D and 3D FE models of the DCB, SLB and ENF specimens are shown in Figure 7 
to Figure 9, which also illustrates the boundary conditions applied for these models. The models 
were based on those presented previously [26]. An additional ENF model with a refined mesh was 
generated, which is shown in Figure 10 with boundary conditions, and discussed further in the 
results section. All models are summarized in Table 2, which lists the analysis solver and dimension 
(2D or 3D) for each model, with a reference to the corresponding mesh. 

Along the specimen length, all models were divided into various sections with different 
mesh refinement. Although the dimensions of the elements varied, all meshes were based on a 
refined mesh region that used elements of 0.5 mm length in the crack growth direction in the region 
immediately around the delamination front. This element length was selected in previous studies 
[28, 31], in which mesh convergence investigations were performed. All models used a uniform 
mesh across the width direction. The specimens with unidirectional laminates used six elements 
through the specimen thickness (2h), while the multi-directional laminate specimens used a refined 
mesh based on a ply level mesh refinement around the delamination interface, as shown in Figure 8. 
Elements away from the delamination interface used a smeared material property for each 
sublaminate based on the rule of mixtures. Further detail and justification of the modeling is given 
in Reference [26]. 
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For all models, the plane of delamination was modeled as a discrete discontinuity in the 
center of the specimen. To create the discrete discontinuity, each model was created from separate 
meshes for the upper and lower part of the specimens with identical nodal point coordinates in the 
plane of delamination. Two surfaces (top and bottom) were created on the meshes as illustrated in 
Figure 6, and a bonded (glued) contact specified between them. For the DCB and SLB specimens, 
the pre-crack was not included in the contact definition as no interaction between the two surfaces in 
this region was expected. For the ENF specimens, contact was required across the pre-cracked 
interface, so the entire surface was included in the bonded contact, and the pre-cracked region was 
created by deactivating the bonded contact for some nodes. For models that were analyzed without 
crack propagation (calculation of strain energy release rate only) as part of developing benchmark 
results as described later, the crack propagation capability was turned off. 

The delaminated surfaces of the DCB and SLB specimens do not come into contact during 
their respective deformation, and thus unbonded contact was not necessary when analyzing these 
specimens. Avoiding the use of contact modeling is advantageous, as the inclusion of a contact 
algorithm significantly increases the solution time of an analysis. In these cases, only bonded 
contact was required in order to simulate delamination growth. The ENF specimen on the other 
hand does involve contact of the delaminated surfaces during deformation, and so unbonded contact 
was necessary to ensure that these surfaces did not interpenetrate during an analysis. 

All specimens were analyzed with 3D models. Additionally 2D models of the DCB and 
ENF specimens were created. A summary of all models run in the current work is presented in 
Table 2. The 3D models used 8-node reduced integration solid brick elements (Marc element type 
1 1 7 [23] and MD Nastran™ element type CHEXA [22]). The 2D models used two-dimensional 4- 
node plane strain elements (Marc™ element type 1 1 and MD Nastran™ element type CQUAD4). 
The 2D and 3D models for a specimen used the same mesh scheme, as shown in Figure 7 and 
Figure 9. The nonlinear solvers in Marc and MD Nastran were used, which applied a full 
Newton-Raphson solution procedure with a load residual tolerance of 0.001. The specimen loading 
was defined in terms of applied displacements, to minimize problems with numerical stability of the 
analysis caused by the unstable propagation. All models were run on a 32-bit Intel Core 2 Duo 2.25 
GHz CPU processor. 

Comparing codes, the input requirements for the VCCT analysis were the same, except for 
the lack of mixed-mode parameters in SOL 400 and SOL 600 as previously discussed. Details on 
the specific keywords and parameters required for a VCCT analysis in each code are given in the 
Appendix. The differences between the codes in terms of the more practical aspects of the 
application are detailed in the Discussion section. 

5. ANALYSIS 

5.1 Computation of Strain Energy Release Rates 

For each of the specimens, the computed strain energy release rate distribution across the 
delamination front was plotted versus the normalized specimen width, y/B, as shown in Figure 1 1 to 
Figure 13, for the 3D and 2D models. Distributions were calculated for a range of different 
delamination lengths, as part of the numerical benchmark results generation described in the next 
section. As strain energy release rate distributions for different delamination lengths were 
previously presented for the DCB and SLB specimens [26], the results in Figure 1 1 and Figure 12, 
respectively, are given for only one delamination length. The ENF results in Figure 13 are given for 
all delamination lengths investigated, for the same applied displacement. For the DCB and ENF 
specimens, only the dominant mode component is shown (I and II respectively), as the others were 
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negligible, and there was no 2D model for the SLB specimen. For the DCB and SLB specimens, 
Figure 11 and Figure 12 demonstrate the excellent agreement between the results in Marc and 
those calculated using Abaqus® in previous work [26], The results for all specimens demonstrate 
that due to 3D effects, such as anticlastic bending of the loading specimen arms, the strain energy 
release rate distribution was non-constant across the delamination front, even for single-mode 
dominated specimens. These curved distributions can cause an initial straight front to grow into a 
curved front. This process has been observed both analytically and experimentally [26, 35-39], 

For the ENF specimen, the results in Figure 13 show that the strain energy release rates 
initially increased with increasing delamination length, which would cause unstable delamination 
growth. This trend changed with respect to the location of the central loading pin, which 
corresponds to a delamination length of 70 mm. The average strain energy release rate initially 
decreased with an increase in delamination length, up to a maximum length corresponding to the 
loading pin location, and then began to increase with longer delamination lengths. This suggests that 
crack growth at delamination lengths past the loading pin was stable. From the results in Figure 13, 
it was found that the distribution of Gu across the delamination front changed slightly with 
delamination length. In this case, Gu peaks on the edges of the delamination front at shorter 
delamination lengths, and peaks in the center at longer delamination lengths. The distribution also 
becomes more curved at delamination lengths approaching the loading pin. The unstable nature of 
delamination growth in ENF specimens loaded quasi-statically, under displacement control, is well 
documented in the literature [40-42], 

5.2 Creating Benchmark Delamination Growth Results 

The approach developed previously for generating numerical benchmark delamination 
growth results [26] was applied to the DCB-1, SLB and ENF specimens. For these specimens, a 
benchmark result set was extracted from a series of models with different delamination lengths. 
These models did not simulate delamination propagation, and were only used to get the load- 
displacement response and the strain energy release rate distribution for different delamination 
lengths. For each delamination length modeled, a failure index was calculated across the 
delamination front using Equation 1, with the Reeder mixed-mode criterion used to compute G c 
(Equation 3). Delamination growth onset was assumed when the failure index at the center of the 
specimen (y/B = 0) reached a value of unity. Specimen load-displacement response up to this point 
was also assumed to be linear. Subsequently, the displacement and load at delamination growth 
onset, dcrtt and P cr n, respectively, was computed by linearly scaling the prescribed displacement and 
load in each analysis (<5 and P respectively), using the following relations [26]: 



The benchmark result set is constructed by plotting the displacement at delamination growth onset 
versus delamination length, as illustrated in Figure 14 for the DCB-1 specimen. This form is used 
for the benchmark result because all specimens were loaded via a prescribed displacement, and the 
change in delamination length is the most appropriate output for assessment of delamination 
growth. In previous work [26], the benchmark results were presented as load-displacement curves, 
to illustrate the application of Equation 5. 
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The benchmark delamination growth results for the DCB-1, SLB and ENF specimens are 
presented in Figure 14, Figure 15 and Figure 16, respectively. Results from previous analyses with 
Abaqus® [26] are also included for the DCB-1 and SLB specimens. The DCB-1 results show that 
the applied opening displacement increased with increasing delamination length, which indicates 
stable growth under displacement control. In contrast, both the SLB and ENF results indicate initial 
unstable delamination growth, as the critical displacement decreases with increasing delamination 
length. The comparison between the Abaqus® and Marc™ results in Figure 14 and Figure 15 shows 
that the two solvers gave almost identical results for these models. 

5.3 Delamination Propagation Analysis: DCB-1 

For the DCB-1 specimen, the previous results in Figure 14 show that delamination 
propagation was predicted to initiate at an applied opening displacement (<5/2) of 0.75 mm. Based on 
this, a two-step loading procedure was applied for the delamination propagation analysis, which 
involved using coarse time increments until just before failure, and fine increments for the region 
involving delamination propagation. Dividing the first step into relatively coarse time increments 
was possible as the load-displacement behavior of the specimen up to failure was expected to be 
linear. In the first step, a prescribed opening displacement of 8/2 = 0.7 mm was applied in 10 
increments. In the second step, the total prescribed displacement was increased to 812 = 1 .0 mm. 
This was applied with a fixed time increment scheme of 50 increments (812 = 0.006 mm each 
increment). 

In both Marc and MD Nastran , the solver has the capability to cut back the increment 
size in the event of convergence issues. In addition, the Marc solver has the capability to activate 
damping when the time step is reduced below a defined minimum. Critically, no convergence issues 
were seen throughout any analyses, and damping was not required. This is quite different behavior 
from that seen previously with Abaqus® [26], where convergence issues associated with 
delamination growth caused significant cutbacks in the time increment, required investigation of 
suitable damping parameters and involved an increase in computational expense. As a result, run 
times for delamination propagation analyses were within a minute for the 2D models and generally 
within a few hours for the 3D models, depending on the selection of increment size and amount of 
delamination growth. The efficiency of the solver is discussed further in a later section. 

The DCB-1 specimen was analyzed with Marc™, MD Nastran™ SOL 400 and MD Nastran™ 
SOL 600, with both 2D and 3D models, as shown in Table 2. All solvers detected delamination 
propagation using simple single-mode criteria, as the SOL 400 and SOL 600 solvers did not have 
any mixed-mode criteria as previously discussed. The results of all analyses are shown in Figure 17 
to Figure 1 9, where Figure 1 7 to Figure 1 8 are the results from only the Marc analysis, and Figure 
19 is a comparison of the delamination growth results from all solvers. 

From the delamination growth results in Figure 17, the VCCT technology gave very close 
comparison with the benchmark results, for both 2D and 3D models. The 3D models showed 
delamination growth at slightly lower applied displacements, which was due to the slightly higher 
strain energy release rates in the 3D models as shown in Figure 1 1 . 

From Figure 18, the delamination was seen to propagate as a straight crack front, which 
contradicts the curved strain energy release rate distribution shown in Figure 11. This propagation 
of a straight front is related to the coarse element size and was also seen in previous analyses [26, 
43], 

From Figure 1 9, all solvers produced almost identical results, for both 2D and 3D models. 
This indicates that the implementation of the VCCT technology is common across solvers. The only 
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discrepancy of note was seen for the 3D model, where on two occasions the SOL 400 model 
showed delamination growth an increment earlier than the Marc and SOL 600 models. This can be 
seen in Figure 1 9 for the first and last delamination growth increments shown (at Aa = 0 mm and 
A a = 4.5 mm). This was attributed to slight numerical differences between the solvers. 

The characteristic step pattern seen in the delamination growth results in Figure 17 and 
Figure 19 was caused by the step change in length as the delamination was grown one element 
length at a time. In the case of the 3D models, the delamination was seen to propagate across the 
crack front within the same increment. This step pattern is illustrated in Figure 20, where the 
numbered stages show the initial increase in applied displacement (1-2), sudden nodal release (2-3), 
increase in displacement until next delamination growth (3-4), and repetition (4-5, 5-6, etc.). These 
step changes in delamination length at a fixed applied displacement would produce a sequence of 
drops in the reaction load and the corresponding “saw-tooth” pattern in the load-displacement 
response [43, 44], The results also show that the comers of the steps, just before a delamination 
growth event, provide the most suitable comparison with the benchmark results, which also allows 
coarse meshes to be adequately used. 

5.4 Delamination Propagation Analysis: DCB-2 

For the DCB-2 specimen [27], the analysis used the model shown in Figure 7, with the 
width across the delamination front increased from 25 mm to 30 mm. The analysis used two load 
steps with coarse and fine time step incrementation, as discussed for the DCB-1 specimen. In the 
first step, a prescribed displacement of <5/2 = 0.9 mm was applied in 2 increments, while in the 
second step the total prescribed displacement was increased to <5/2 = 6.0 mm in 250 increments 
(<5/2 = 0.0204 mm each increment). As with the DCB-1 specimen, no significant convergence 
issues were seen and damping was not required. 

The DCB-2 specimen was analyzed in Marc with both 2D and 3D models as shown in 
Table 2. Delamination propagation was detected using single-mode criteria. Numerical results were 
also taken from a previous publication [45], in which the DCB-2 specimen was analyzed using 
Abaqus® 6.8. Further detail on the modeling approach adopted in that work is given in References 
[21, 26, 45], which includes discussion on the damping and other parameters required. Benchmark 
analytical solutions for delamination growth were taken from Reference [27], where relationships 
between the applied displacement and delamination length are presented based on beam theory. 
Experimental delamination growth results were taken from Reference [46], where the delamination 
length was measured at intervals on the edge of the specimen. The experimental results were 
included since the original experiments showed very little fiber bridging [46], which is not 
accounted for in the analysis. The delamination growth results of all analyses and benchmarks are 
shown in Figure 21 . 

It is observed from Figure 21 that the VCCT technology was capable of accurately 
representing the analytical and experimental benchmark results in terms of delamination growth 
behavior. Excellent agreement was seen between the two analysis codes, and between the numerical 
analyses and both the experimental and analytical results. The change in character of the Marc 
results at a delamination growth Aa = 40 mm was due to the delamination reaching a region of the 
mesh with larger elements, which led to larger step changes in delamination growth. At the 
transition between the fine and coarse mesh regions, the Marc results showed a deviation from 
the expected behavior. As discussed previously for the DCB-1 specimen and illustrated in Figure 
20, the peak displacements just prior to delamination propagation in each step give a consistent 
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alignment with numerical benchmark results. The results at the mesh transition region indicated 
in Figure 21 demonstrate an error in the VCCT implementation, caused by the difference between 
element lengths ahead and behind the delamination front. The Abaqus® results shown in Figure 21 
do not show this error since a consistently fine mesh was used. However, this error was seen in 
other Abaqus® results [45], and indicates a similar implementation error in both codes. The issue of 
mesh transition is further investigated and discussed in detail in Section 6. 1 . 

The results for the DCB-2 specimen demonstrate the difference between code 
benchmarking using experimental and analytical approaches, in comparison with the use of the 
numerical benchmark results for the DCB-1. The use of experimental and analytical results is 
highly valuable in validating that the simulation is capable of capturing the physical phenomena. 
These approaches, however, may introduce additional factors that can influence an assessment of 
simulation capabilities, which include statistical variance, the use of approximations, or 
complications such as fiber bridging. As these aspects were not expected to be significant in this 
work, the comparisons with benchmark results for the DCB-1 and DCB-2 specimens produced 
similar conclusions. Further, the results for the two specimens demonstrate that numerical 
benchmarks are suited for efficiently isolating and assessing code capabilities as they involve a 
direct comparison between the behavior of identical models. The use of numerical benchmarks 
would also be highly valuable where experimental and analytical results are unavailable or 
difficult to obtain, such as for complex specimens and structures. 

5.5 Delamination Propagation Analysis: SLB 

For the SLB model, shown in Figure 8, the analysis used two load steps with coarse and 
fine time step incrementation as discussed for the DCB-1 specimen. In the first step, a prescribed 
center displacement of w = 3.0 mm was applied in 6 increments, while in the second step the 
total prescribed displacement was increased to w = 4.0 mm in 100 increments 
(w = 0.01 mm each increment). As with the DCB specimens, no significant convergence issues 
were seen and damping was not required. 

One aspect that was critical for the analysis of the SLB specimen was that the crack growth 
direction needed to be specified with a user-defined vector, instead of with the default automatic 
algorithm that, as previously discussed, is based on a principal stress criterion. In preliminary SLB 
models using the default approach, errors were found in the propagation of the delamination front, 
where the solver would detect the critical node but not be able to locate the appropriate node in the 
intact region to grow the delamination front. This was found to be due to the calculated crack 
growth direction, which in some instances pointed in directions from which no suitable nodes could 
be found. This problem was also found to be particularly critical for edge nodes, and was found to 
occur for both coarse and fine meshes. In instances where the delamination could no longer be 
propagated due to the inappropriate crack growth direction, the entire delamination was made 
“inactive”, and all VCCT calculations were ceased. 

In order to address this problem, it was necessary to specify a user-defined vector for the 
crack growth direction, which was possible for the SLB specimen, as the crack growth direction 
was known. For the SLB specimen as shown in Figure 8, the crack growth direction was aligned 
with the global x-direction (i.e. the user-defined crack growth direction vector was set to 
<1 ,0,0>). It should be stated that the calculation of the strain energy release rates is not based on the 
crack growth direction, but uses a local crack front coordinate system as previously described. As 
such, setting the crack growth direction to be fixed acts only to assist the solver in locating the 
correct node for crack propagation, and does not influence the mode mixity at the crack front. 
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The consequence of this issue was that SOL 400 and SOL 600 could not be used to analyze 
the SLB specimen, as only the default crack growth direction option was available to be selected in 
the version of MD Nastran appied in this work. More critically, the SLB specimen is a mixed- 
mode problem, so that in any event, as only single-mode crack growth could be implemented, SOL 
400 and SOL 600 could not have been used. 

The SLB specimen was analyzed in Marc with a 3D model as shown in Table 2. 
Delamination propagation was detected using the Reeder mixed-mode criterion. A 2D model was 
not analyzed, as the multi-directional laminate could not be modeled correctly. The delamination 
growth results are shown in Figure 22, and the progression of the delamination front is illustrated in 
Figure 23. 

From these results, the VCCT technology was again capable of accurately representing the 
delamination growth benchmark results for this mixed-mode case as shown in Figure 22. The 
benchmark results show that the SLB specimen exhibited an initial period of unstable delamination 
growth. This was caused by the strain energy release rates remaining critical with increasing 
delamination length. In the benchmark curve, unstable delamination growth is indicated by a 
reduction in applied displacement below the initiation point. As the analysis was run in 
displacement control, it was not possible to get a reduction in displacement, so the unstable 
delamination growth is seen as a large crack growth step. This corresponded to the delamination 
extending by 21 mm within one increment, which is shown in both Figure 22 and Figure 23. This 
was followed by stable delamination growth, where the benchmark and FE results correlated very 
closely. The delamination growth behavior as shown in Figure 23 demonstrated that a curved 
delamination front developed and propagated. This was caused by the curved strain energy release 
rate distribution as shown in Figure 12. As described for the DCB-1 specimen, this led to a slight 
difference compared to the benchmark results, which were generated based on a straight 
delamination front. It is possible that the exact shape of the delamination front shown in Figure 23 
was affected by the delamination propagating through a mesh transition, which is further 
demonstrated for the ENF specimen and investigated in detail in the Discussion section. 

5.6 Delamination Propagation Analysis: ENF 

For the ENF model shown in Figure 9, the analysis used two load steps with coarse and 
fine time step incrementation as discussed for the DCB-1 specimen. In the first step, a prescribed 
center displacement of w = 4.0 mm was applied in 2 increments, while in the second step the 
total prescribed displacement was increased to w = 8.0 mm in 100 increments (w = 0.04 mm each 
increment). As with all previous specimens, no significant convergence issues were seen and 
damping was not required. However, as found for the SLB specimen, the default crack growth 
direction algorithm resulted in delamination propagation errors, and a user-defined vector was 
used. Again, due to the requirement for a user-defined crack growth vector, MD Nastran SOL 
400 and SOL 600 could not be used. 

The ENF specimen was analyzed in Marc with both 2D and 3D models as shown in 
Table 2. Delamination propagation was detected using single-mode criteria. The results of the 
analyses are shown in Figure 24 and Figure 25. As observed from these results, the VCCT 
technology was capable of accurately representing the benchmark delamination growth results for 
the baseline ENF mesh. As with the DCB-1 specimen, delamination growth was initiated in the 3D 
model at a slightly lower applied displacement than the 2D model, which was due to a slightly 
higher strain energy release rate as shown in Figure 13 for a = 30 mm. 
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The results for the ENF specimen showed initial unstable delamination propagation 
followed by stable propagation, in the same way as described for the SLB specimen. The first 
delamination growth event corresponded to a delamination growth step of 41 mm, as shown in 
Figure 24. The unstable delamination growth was followed by stable growth, where the benchmark 
and FE results correlated very closely. As previously discussed, the change from unstable to stable 
delamination growth is associated with the inflection point of the applied displacement versus 
delamination growth behavior, which is indicated in Figure 24. Although in experimental testing 
focused on material characterization, the delamination is not extended past this point, for the 
numerical analysis this region remains valuable for benchmarking purposes. 

From Figure 25, the delamination front formed in the baseline mesh by the unstable 
delamination growth was jagged and non-straight across the width. This irregular and unexpected 
pattern was considered a product of the delamination propagating through a mesh transition region, 
in addition to insufficient mesh density across the width. A refined model, shown in Figure 10, was 
analyzed in which the delamination propagation only occurred in the fine mesh region, where the 
element length was 0.5 mm in the crack growth direction. Rather than using a large fine mesh 
region to capture all of the unstable delamination growth from a = 30 mm to a = 71 mm, it was 
decided to instead increase the length of the initial delamination to 65 mm. This located the 
delamination within the unstable delamination growth region, so that the stable delaminaiton growth 
behavior would still be identical to the baseline case, while keeping the model size reasonable. 
Additionally, the mesh density across the width was doubled, while for efficiency the length of the 
fine mesh region was reduced. 

The results of the analysis of the ENF refined mesh model are shown in Figure 26, where 
the increased mesh density and initial delamination of 65 mm produced a straight delamination front 
without the jaggedness across the width. The delamination growth results for the models are shown 
in Figure 27, where the initial delamination length of a = 65 mm is represented as a delamination 
growth Aa = 35 mm, so that the baseline for delamination growth in the two numerical analyses is 
the same. From the delamination growth results, it can be seen that despite different delamination 
front patterns, the models gave very close comparison with the benchmark curve for the stable 
delamination propagation stage. These results demonstrate that crack propagation through a mesh 
transition region could lead to an incorrect delamination front developing and propagating, and that 
the mesh refinement level was critical to accurately capturing the shape of the delamination front. 
The results in Figure 27 also clearly demonstrate the effect of large elements on the delamination 
growth, where the peaks of the large steps just prior to delamination growth gave the most suitable 
comparison with the benchmark results. 


6. DISCUSSION 
6.1 Mesh transition 

The results in Figure 21 and Figure 25 demonstrated an error in the implementation of the 
VCCT in Marc , with regards to calculation of strain energy release rates, where elements ahead 
and behind the delamination front had different lengths. This was also seen for the Abaqus® solver 
in Reference [45]. As detailed in References [19, 20], the VCCT equations require a modification 
factor to account for uneven element lengths ahead and behind the crack front. The results in 
Figure 21 and Figure 25 demonstrate that this has not been implemented, so that uneven element 
lengths, such as those seen as a crack propagates through a mesh transition region, result in incorrect 
calculation of the strain energy release rate. To illustrate this effect further, a 2D DCB-1 model with 
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several mesh transition regions was analyzed in Marc™, with the mesh and the results shown in 
Figure 28. The errors are most clear for the transition to smaller elements, as is shown at 
A a = 5.25 mm, where an under-estimation of the strain energy release rates due to the mesh 
transition prevents delamination growth from occurring and following the benchmark results. The 
same error is less evident though still occurring in the two transitions to larger meshes at Aa = 1 .0 
mm and A a = 7.25 mm. For the 2D DCB-1 model in this example, the errors did not propagate, so 
that predictions in a regular mesh following a mesh transition region returned to follow the 
benchmark solution. However, the results for the 3D ENF model demonstrated that an incorrect 
delamination front could develop due to a mesh transition region and be propagated. 

6.2 Solver comparison: Marc™ and Abaqus® 

Comparing Marc results in this work (which includes the adaptation of the Marc 
implementation into MD Nastran™) with previous Abaqus® results [26], one clear difference was 
the solution process. Critically, in the Marc implementation of VCCT [23], once crack growth is 
detected in an increment the crack front node is released and the increment is restarted. This allows 
for multiple crack growth instances to occur in one increment, and allows for coarse time 
increments to be used, which is computationally efficient. By comparison, in the Abaqus® 
implementation [21], only one crack growth instance can occur in each increment, so it is necessary 
to reduce the increment size considerably in order to ensure that no over-estimation occurs. The 
failure index determined from the strain energy release rates is monitored, and the user defines a 
limit on the amount that the failure index can exceed a value of 1 .0. This overshoot limit, or “release 
tolerance”, as such becomes another parameter in the model that requires careful selection [26], So, 
while the Abaqus® approach decreases the increment size to suit the crack growth, in Marc™ the 
crack growth is increased to suit the increment size, so that larger increments can be used with 
increased computational efficiency. 

The convenience of using large increments needs to be managed carefully with the need to 
capture the initiation point for crack growth. This is illustrated in Figure 29, where a DCB-1 
specimen is analyzed in Marc with fine and coarse time incrementation. With fine time 
increments, a step pattern of crack growth is observed, which indicates that no over-estimation is 
occurring due to the increment size. The analysis with coarse time increments is seen to follow the 
fine solution closely, though the use of coarse increments means that the exact initiation point for 
crack growth is not captured. This issue can be addressed by using different time incrementation 
throughout the analysis, where for example coarse increments are used until just before crack 
growth initiation (which could be determined from a preliminary coarse analysis), fine increments 
are used to capture the crack growth initiation, and if necessary coarse increments can be re-used to 
capture crack propagation. This is similar to the two-step approach applied in this work. 
Alternatively, it may be considered that the benefit of using coarse increments outweighs the slight 
increase in accuracy in capturing the initiation point, which may be more appropriate for analysis of 
large structures. 

Another important difference between Marc™ and Abaqus® is the convergence difficulties 
and the subsequent damping required. In previous work using Abaqus® [26], it was found that 
damping, or “stabilization”, needed to be added to the solver, in order to get a solution in light of the 
convergence issues. The introduction of stabilization parameters requires considerable effort in 
parametric investigation in order to determine a suitable compromise between damping and solution 
accuracy. In contrast, for the models considered in this work with Marc , no severe convergence 
issues were recorded and damping was not applied in any of the solutions. 
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The two aspects of differing time increment requirements (caused by allowing multiple 
crack growth events in an increment) and the requirement for damping meant that the Marc solver 
was considerably more computationally efficient than the Abaqus® solver for the models analyzed. 
This was most noticeable for the 3D models, where for example, the Abaqus® solver required run 
times of up to several days [26], while the Marc solver experienced run times of only several 
hours. Although the analyses were run on different machines, the comparison clearly highlights the 
considerable difference between the two solvers. This difference is further exacerbated by die 
introduction of the release tolerance and stabilization parameters in Abaqus®, which typically 
require parametric investigations to determine. 

One aspect that required attention in the Marc analysis was the definition of the crack 
growth direction, where the default automatic algorithm was unable to propagate the crack 
correctly. With the user-defined vector, it was possible to achieve the correct solution for the 
problems considered. However, for problems where the crack growth direction is unknown or 
changes along the delamination front, this may present further difficulties. While information on the 
exact procedure implemented in Abaqus® is not available in the documentation, none of the 
previous analyses of DCB and SLB specimens indicated any issues with locating the correct nodes 
for crack propagation [26, 45]. 

The use of the crack growth direction in Marc is an important difference from the 
implementation in Abaqus®, and has both advantages and disadvantages. In Abaqus®, the crack 
front node is simply “released” from contact, and the new crack front is found by determining 
which nodes lay on the boundary of the contact definition. This is a simple technique that ensures 
that there are no issues with propagation in terms of locating the correct nodes, which are located 
implicitly. In Marc , the use of the crack growth direction as a search vector is disadvantageous, as 
it introduces the possibility that either incorrect nodes are found, or no nodes are found at all. The 
difficulty in determining the correct crack growth direction is likely to be caused by the application 
of the principal stress criterion, where the accuracy of this criterion at the interface between two 
different orthotropic materials is uncertain, as the method was originally developed for isotropic 
metals. 

On the other hand, the use of the crack growth direction in Marc has the advantage of 
allowing the crack growth to be explicitly tracked for each node. As the crack front moves from one 
node to another, information regarding the accumulated crack growth is directly passed between 
nodes. In this way, each crack front node at any increment maintains the total crack growth it has 
taken for the crack to reach that node since the beginning of the analysis. This makes it very easy to 
track the crack growth at each node across the crack front, and also to implement R-curve behavior 
where the fracture toughness parameters vary with crack length. 

Despite the advantages in explicitly tracking the crack growth data, it is possible to achieve 
a determination of the total crack growth for any node and a capacity to implement R-curve 
behavior without requiring the use of the crack growth direction. For R-curve behavior, this can be 
achieved by defining a spatially varying field, where the fracture parameters are set as a function of 
geometric coordinates. In the examples in this work, the fracture parameters could be made 
dependent on one coordinate (the x-coordinate for all specimens as shown in Figure 3 to Figure 5) 
though this would change depending on the configuration of the specimen and the pre-crack. This 
approach would be possible within Abaqus®, where fracture parameters can be made dependent on 
a state variable, which can vary according to a given field. Separately, for tracking crack growth 
data at each node, this can be achieved by determining the minimum distance between a crack front 
node and any node defined as part of the initial crack front. This capability is currently not 
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implemented into either code, though could easily be done during post-processing or with a user 
subroutine. 

More broadly, there are other differences between the implementation of the VCCT into 
Marc™ and Abaqus®. Both codes implement the crack propagation modeling as a disconnection of a 
contact surface between two initially bonded or “glued” interfaces. However, while in Abaqus® this 
implementation is only available for use with a contact definition, in Marc it is possible to apply 
this to multi-point constraints. Furthermore, in Marc , the VCCT can be used to model crack 
propagation by remeshing or by disconnecting elements along their edges, though at this stage these 
techniques are only available for 2D solid (plane stress/strain) or shells. Marc also offers a user 
subroutine “ucrackgrow.f ’, so that the crack growth behavior using the VCCT can be customized to 
specific requirements. This allows the user to control the crack growth direction, which is used to 
locate where the crack will propagate, and the amount of crack growth in an increment. 

6.3 Solver comparison: Marc™ and MD Nastran™ 

Comparing Marc™, MD Nastran™ SOL 400 and SOL 600, the results presented in this work 
indicate that the implementation of the VCCT is identical. However, as discussed previously, in 
general the range of technologies and options implemented into MD Nastran tends to be limited, 
as it can take some time and several releases for aspects developed natively in Marc to be migrated 
within the MD Nastran environment [34], As such, comparison of the codes is simplified to 
determining which of the technologies and options that are in Marc are actually available within 
MD Nastran , and understanding what additional options need to be added to an MD Nastran 
model to apply the VCCT technology. Critically, in this work MD Nastran could not be used for 
any mixed-mode delamination growth specimens such as the SLB specimen as there were no 
mixed-mode criteria implemented. Additionally, SOL 400 and SOL 600 could not be applied with 
the SLB and ENF specimens due to the difficulties in using the default algorithm for determining 
the crack growth direction. 

In terms of SOL 400, the VCCT implementation is integrated within the MD Nastran™ 
nonlinear solver. In the analysis in this work, and in previous benchmark studies [47], the MD 
Nastran nonlinear solver was more robust and efficient than in Marc This is partly due to a 
difference in the control of the increment size. In MD Nastran , this is based on the convergence 
parameters, while in Marc control is framed around either the number of iterations in an increment 
or the damping parameter. Additionally, the SOL 400 solver in the version of MD Nastran applied 
in this work (v2008.0) could not be used with user subroutines in order to customize the solver 
behavior. In Marc , the “ucrackgrow.f’ subroutine can be used to define the crack growth direction 
and amount of crack growth, as previously discussed, and there are a wide range of other 
subroutines available to customize different aspects of the solver [23]. It was also seen that the SOL 
400 had larger memory requirements, and generated larger output files and larger temporary files 
while running. 

For SOL 600, as the solution runs natively in Marc , care needed to be taken with the MD 
Nastran input deck to ensure that the model was translated into Marc as intended. This required 
ensuring aspects such as the material orientations, element types and output requests were translated 
correctly, and that the cards controlling the VCCT technology were included. This was particularly 
challenging where aspects of one code did not easily “map” or translate into the other, or where 
access to Marc functionality was limited by the keywords available in the MD Nastran 
framework. To this end, an understanding of the Marc input file and solver technologies was 
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considerably beneficial in ensuring that the MD Nastran™ input file was translated into a Marc™ 
input file with the appropriate format and content. 


7. SUMMARY AND CONCLUSIONS 

The delamination propagation simulation capabilities of the commercial FE codes Marc 
2008rl and MD Nastran™ v2008.0 with the VCCT were assessed. Benchmark delamination growth 
results for DCB, SLB and ENF specimens were generated using a previously developed numerical 
approach [26]. Experimental and analytical benchmark delamination growth results were taken 
from the literature for a second DCB specimen. Specimens were analyzed using 3D and 2D models, 
and compared to previous analyses using Abaqus® with the VCCT implemented. The results 
showed the following: 

• The approach previously developed for using numerical benchmarks is a practical and 
efficient method for assessing delamination growth behavior. This approach offers 
advantages over benchmarking using experimental and analytical results for isolating 
code capabilities as it involves a direct comparison between identical models. 

• The VCCT implementation in Marc and MD Nastran was capable of replicating the 
benchmark data for the DCB, ENF and SLB specimens, in terms of delamination growth 
behavior. 

• For the Marc™ and MD Nastran® analyses, no issues were seen with convergence, so 
that no damping was applied, and in general, the analyses were not considered highly 
sensitive to the solver input parameters. 

• The capacity in Marc and MD Nastran to automatically internally restart the 
increment after crack growth, and hence capture multiple crack growth events in one 
increment, was highly beneficial in terms of computational efficiency. 

• Based on the previous two points, the analyses in Marc and MD Nastran were 
considerably more computationally efficient than those previously conducted in 
Abaqus®. Problems involving identical meshes across all codes took several hours in 
Marc™ and MD Nastran™, and several days in Abaqus®, although these were run on 
different machines. 

• The implementation of VCCT in Marc™ and MD Nastran® uses a principal stress 
criterion to automatically determine the crack growth direction and hence explicitly 
locate the nodes required for crack propagation. This procedure was problematic for the 
problems considered and was not required in the Abaqus® analyses. 

• For the SLB and ENF specimens, problems with the default automatic algorithm for 
determining crack growth direction in Marc™ and MD Nastran® were overcome by 
manually defining a suitable vector. 

• The implementation of the VCCT in MD Nastran is a nominally identical, yet limited 
version of that available in Marc , where the absence of mixed-mode crack growth and 
the option to specify a user-defined crack growth direction were critical in this work. 

• The implementation of the VCCT in Marc and MD Nastran does not account for 
irregular mesh lengths ahead and behind the crack front. 

• The use of mesh transition and fine mesh regions needed to be managed carefully in 
order to ensure an appropriate delamination front was formed. 
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Overall, it is clear that delamination propagation modeling is rapidly evolving in 
commercial FE codes, with each new software release introducing new technologies as well as 
refinements and extensions on previous technologies. This highlights the need for benchmarking 
techniques and examples that are capable of isolating and assessing the key requirements for 
delamination propagation simulation. 
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TABLE 1 . MATERIAL PROPERTIES . 


T300/1076 Unidirectional Graphite/Epoxy Prepreg 

En = 139.4 GPa 

E 22 = 10.16 GPa 

E 33 = 10.16 GPa 

vi 2 - 0.30 

vi 3 = 0.30 

V23 - 0.436 

G 12 = 4.6 GPa 

G 13 = 4.6 GPa 

G 23 = 3.54 GPa 

C12K/R6376 Unidirectional Graphite/Epoxy Prepreg 

En = 146.9 GPa 

E 22 = 10.6 GPa 

E 33 = 10.6 GPa 

vi 2 = 0.33 

vi 3 = 0.33 

V23 - 0.33 

G 12 = 5.45 GPa 

G 13 = 5.45 GPa 

G 23 = 3.99 GPa 

T800/924 Unidirectional Graphite/Epoxy Prepreg 

Ei i = 126 GPa 

E 22 = 7.5 GPa 

E 33 = 7.5 GPa 

v\ 2 = 0.263 

vi 3 = 0.263 

V23 - 0.263 

Gi 2 = 4.981 GPa 

G 13 = 4.981 GPa 
G/ c = 0.281 kJ/m 2 

G 23 = 3.321 GPa 


The material properties are given with reference to the ply coordinate axes where index 1 1 denotes the ply 
principal axis that coincides with the direction of maximum in-plane Young’s modulus (fiber direction). Index 22 


denotes the direction transverse to the fiber in the plane of the lamina and index 
perpendicular to the plane of the lamina. 

TABLE 2. FINITE ELEMENT MODELS. 

33 denotes the direction 

Model 

FE Solver 

Dimension 

Reference 

Benchmark 

DCB-1 

Marc™ 

2D 

Figure 7 a 

Numerical - Marc™ 



3D 

Figure 7b 


MD Nastran™ SOL 400 

2D 

Figure 7 a 




3D 

Figure 7b 


MD Nastran™ SOL 600 

2D 

Figure 7a 




3D 

Figure 7b 



Abaqus® 

3D 

Figure 7b [26] 

Numerical - Abaqus® 

DCB-2 

Marc™ 

2D 

Figure 7a 

Experimental [46] and 



3D 

Figure 7b * 

analytical [27] 


Abaqus® 

2D 

Figure 7a [45] 




3D 

Figure 7b [45] 


SLB 

Marc™ 

3D 

Figure 8 

Numerical - Marc™ 

ENF 

Marc™ 

2D 

Figure 9a 

Numerical - Marc™ 



3D 

Figure % 


ENF refined mesh 

Marc™ 

3D 

Figure 10 

Numerical - Marc™ 


* Width increased from model shown, from 25 mm to 30 mm, by duplicating elements, so that the mesh length 
in the width direction remained unchanged. 
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APPENDIX 

Example input for VCCT analysis in Marc™ and MD Nastran™ 

The input required to perform VCCT analysis in Marc™, MD Nastran™ SOL 400 and SOL 
600 is discussed in the following paragraphs. This is focused exclusively on the keywords and 
parameters of relevance for VCCT analysis, where for further detail reference should be made to the 
software documentation [23, 22], Note that the ellipses are used to indicate removed text. Also, 
in most cases example values are included as taken from an input file, and these should not be 
implied as fixed, recommended or necessary unless otherwise stated. 

Example input for Marc™ 

The key parameters required for VCCT analysis in Marc are listed below. Note that in this 
work input files were completely generated using the Mentat pre-processor. 

Model definition options 

Crack front node set. This is useful for the VCCT card. 

define ndsq set crackl_nodes 

8572 8580 8588 8596 8644 


Contact body. The properties that control the contact and the two bounding bodies (groups of 
elements) that form the crack interface. 

contact 

2,0 

• • • 

$ . . . . contact body 1 : cbody_top 

• • • 

$ . . . . contact body 2 : cbody_bot 


VCCT parameters : Parameters that control the VCCT implementation, 
vcct 

1 

1 2 

crackl 

crackl_nodes 
2, <a>, 0, <b> 

0.0, 0.0, <GIc>, <GIIc>, <GIIIc>, <nl>, 0.0, 0.0, ... 

• • • 

<dl> , <d2> , <d3> 
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where <a> is the code for crack growth direction method (1 for principal stress criterion, 4 for 
a specified vector), <b> is the code for crack growth criterion (2 for single-mode, 4 for 
Reeder), <GIc>, <GIIc>, <GIIIc> are the fracture toughness values, <nl> is the Reeder 
exponent used if <b>=4, and <dl>, <d2>, <d3> are the jc, y, z components of the crack 
growth direction if <a>=4. 

Note that it is not required to tell the code which contact interfaces are associated with the 
crack. The user only inputs the crack front nodes, and the solver will determine which contact 
interface these nodes are associated with. 

History definition options 

Multiple load cases'. Useful to perform a multiple-step analysis to save computation time when the 
initiation point for crack growth is approximately known (use the standard loadcase card). 

Nodal output: Useful output to plot in order to visualize the crack front data, where the crack front 
images in this report were generated primarily using the “glue deactivation status”. 


post 




• • • 

nodal. 

38 


contact status 

nodal. 

58 

-» 

GI 

nodal. 

59 


GII 

nodal. 

60 


GUI 

nodal. 

73 

-» 

glue deactivation status 

nodal. 

74 

-» 

VCCT failure index 


Example input for MD Nastran™ SOL 400 and SOL 600 

Examples of the parameters required for conducting VCCT analysis in MD Nastran SOL 
400 and SOL 600 are given below, which follow the comments given previously. In this work, the 
input decks were generated using MSC Patran , with minor editing in order to add parameters 
necessary for VCCT. These are specified below, which also includes other parameters of relevance 
to VCCT. In general, the input files for SOL 400 and SOL 600 are mostly identical, except for a few 
keywords. The comments to the right of the page indicate which parameters are exclusive for each 
solver, and where these are absent the parameters relate to both solvers. 

Executive control statement 


SOL 400 

SOL 600, NLSTATIC OUTR=xdb 


— > SOL 400 only 
— »> SOL 600 only 
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Case control 


VCCT = 1 


—* reference to bulk data entry 


Multiple load cases : Useful to perform a multiple-step analysis to save computation time when the 
initiation point for crack growth is approximately known (use the standard SUBCASE card). 


SUBCASE 1 
STEP 1 

ANALYSIS = NLSTATIC 


SOL 400 only 
SOL 400 only 


SUBCASE 2 
STEP 2 

ANALYSIS = NLSTATIC 


SOL 400 only 
SOL 400 only 


BEGIN BULK 


Contact body. The parameters relating to the glued contact definition. 
BCPARA 0 NLGLUE 1 NBODIES 2 


BCTABLE 0 1 

SLAVE 2 0 . 0 , 

2 0 0 

FBSH 1.+20 .9 

MASTERS 1 


0 , 

0 , 


BCTABLE 1 1 

SLAVE 2 0 . 0 , 

2 0 0 

FBSH 1.+20 .9 

MASTERS 1 


0 , 

0 , 


VCCT parameters'. Parameters that control the VCCT implementation. 
VCCT, 1, 1, 2 

, <GIC> , , , , , , <GIIC> , <GIIIC> 

• • • 

<nodes> 
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where <GIc>, <GIIc>, <GIIIc> are the fracture toughness values, and<nodes> is either 
a sequence of node numbers or the negative index of a node set (to avoid confusion with node 
numbers). 


Crack front node set. Example of a node set, here assigned number “9”, which would be referenced 
as “-9” in the VCCT statement above. 

SET3 9 GRID 8572 8580 8588 


Nonlinear solution parameters : In SOL 400 and SOL 600, the NLPARM card is used. In this work, 
in addition to the NLPARM card, for SOL 400 the NLADAPT card was used, and for SOL 600 the 
NLAUTO and NLSTRAT cards were used. This was done in order to ensure that the nonlinear 
solution scheme was consistent between solvers as much as possible. 


Element properties : In SOL 400 only, an additional element property card is required to access 
advanced properties required for VCCT. Examples are given below for 3D solid elements (which 
reference the PSOLID card) and 2D plane strain elements (which reference the PLPLANE card). 


PSOLID 1 

PSLDN1 1 

C8 

PLPLANE 1 

PSHLN2 1 

C4 


1 1 

1 

SOLID L 

1 

1 1 

PLSTRN L 


REDUCED 


25. 


Contact body. The two bounding bodies (groups of elements) that form the crack interface (here an 
example for 3D, for 2D replace “2D” with “3D” in the third parameter of BCBODY). 

$ Deform Body Contact LBC set: cbody_top 
BCBODY 1 3D DEFORM 1 0 

BSURF 1 385 386 387 

• • • 

$ Deform Body Contact LBC set: cbody_btm 
BCBODY 2 3D DEFORM 2 0 

BSURF 2 385 386 387 
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Interlaminar sliding shear 
Mode II 



Interlaminar scissoring shear 
Mode III 


Figure 1 . Fracture modes. 
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800 


G , 

c 

J/m 2 


G , 

c 

J/m 2 



0 

experimental data 

• 

mean values 


curve fit: G c = G /c+ ^-G^G/G/ 
r/= 1.62 


—I I I— 


_l 1 1_ 


-I I I— 


-I I I— 


-I I I— 


0.2 0.4 0.6 

mixed mode ratio G // /G r 
a. T300/914C 


0.8 



mixed mode ratio GJG T 
b. C12K/R6376 


Figure 2. Mixed-mode fracture criterion for composite materials. 
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A* 



DCB-2 

30.0 mm 

3.0 mm 

150.0 mm 

30.0 mm 

T800/924 

unidirectional 



6 25.4 mm 

2.032 mm 
t 2 2.032 mm 

2 L 177.8 mm 

a 34.29 mm 

C12K/R6376 


[±30/0/-30/0/30/0 4 /30/0/-30/0/-30/30 // -30/30/0/30/0/ -30/0 4 /-30/0/30/0/±30] 

Figure 4. SLB specimen. 



B 25.0 mm 

2 h 3.0 mm 

L 70.0 mm 

2 L 140.0 mm 

a 30.0 mm 

T300/1 076 
unidirectional 
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initial 

delamination 


intact 

region 



a. Specimen definition 


bonded contact 
contact surfaces 


/ 




unbonded bonded contact 

contact contact surfaces 



b. Contact in bonded region only 


c. Contact across entire interface with 
bonded and unbonded contact regions 


Figure 6. Contact-based crack modeling. 



b. 3D model 


Figure 7. DCB-1 deformed specimen mesh with boundary conditions. 
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detail of local 
model around 
delamination front 


[ 30 / 0 /- 30 / 0 /± 30 ] 
[0] 4 

[ 0 / 30 / 0 / - 30 ] 
[ 30 ] 
[- 30 ] 



bonded 

contact 

interface 


[± 30 / 0 /- 30 / 0 / 30 ] 
[0] 4 

[ 30 / 0 /- 30 / 0 ] 


[± 30 / 0 /- 30 / 0 / 30 / 0 4 / 30 / 0 /- 30 / 0 /- 30/30 // - 30 / 30 / 0 / 30 / 0 / - 30 / 0 4 /- 30 / 0 / 30 / 0 /± 30 ] 


Figure 8. Deformed mesh of SLB specimen and detail of region around delamination front with 

boundary conditions. 
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unbonded bonded contact 

contact interface interface 



x, u 


a. 2D model, deformed 



b. 3D model, undeformed 

Figure 9. ENF specimen mesh with boundary conditions. 
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y/B 

Figure 1 1 . Strain energy release rate distribution, DCB (a - 30.5 mm, d/2-1.0 mm). 
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— — Abaqus Benchmark Gl [26] 

o 

Marc 3D Gl 

Abaqus Benchmark Gil [26] 

□ 

Marc 3D Gil 

— — - Abaqus Benchmark Gill [26] 

A 

Marc 3D Gill 



- 0.5 - 0.4 - 0.3 - 0.2 - 0.1 0.0 0.1 0.2 0.3 0.4 0.5 

y/B 

Figure 12. Strain energy release rate distribution, SLB (a - 31 mm, w - 2.8 mm). 
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0.6 
0.5 
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kJ/m 2 

0.3 
0.2 
0.1 
0.0 

- 0.5 - 0.4 - 0.3 - 0.2 - 0.1 0.0 0.1 0.2 0.3 0.4 0.5 

y/B 

Figure 13. Strain energy release rate distribution, ENF (w - 4.0 mm). 
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applied opening displacement 8/2, mm applied opening displacement &2, mm 



delamination growth A a, mm 

Figure 14. DCB-1, benchmark delamination growth results. 



delamination growth A a, mm 

Figure 15. SLB, benchmark delamination growth results. 
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applied center displacement w, mm 


9.0 



0 10 20 30 40 50 60 

delamination growth A a, mm 

Figure 16. ENF, Marc benchmark delamination growth results. 
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Figure 17. DCB-1 , applied displacement versus delamination growth. 
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applied opening displacement d/2, mm 


bonded contact 
interface at 
d/2 = 1 .0 mm 


delamination front 
at 8/2 = 1 .0 mm 



initial delamination 
front location 

Figure 18. DCB-1 bottom sublaminate deformed mesh at <5/2 -1.0 mm, with delamination 

interface. 



Figure 19. DCB-1, applied displacement versus delamination growth. 
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applied opening displacement & 2 , mm 


a 






Figure 20. Delamination growth behavior. 



delamination growth A a, mm 

Figure 21 . DCB-2, applied displacement versus delamination growth. 


38 


applied center displacement w, mm 



Figure 22. SLB, applied displacement versus delamination growth. 



Figure 23. SLB bottom sublaminate mesh with delamination front progression. 
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applied center displacement w, mm 



delamination growth A a, mm 

Figure 24. ENF, applied displacement versus delamination growth. 



Figure 25. ENF bottom sublaminate mesh with delamination front progression. 


40 


applied center displacement w, mm 



Figure 26. ENF refined model bottom sublaminate mesh with delamination front progression. 



delamination growth A a, mm 

Figure 27. ENF, applied displacement versus delamination length. 
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Figure 28. Effect of mesh transition regions, DCB-1 2D. 
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using Abaqus®. The results demonstrated that the VCCT implementation in Marc™ and MD Nastran™ was capable of accurately 
replicating the benchmark delamination growth results and that the use of the numerical benchmarks offers advantages over benchmarking 
using experimental and analytical results. 
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